import pandas as pd
import numpy as np
from scipy.stats import norm
import os
import matplotlib.pyplot as plt
# Black-Scholes Function
def black_scholes(S, K, T, r, sigma, option_type="call"):
d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
d2 = d1 - sigma * np.sqrt(T)
if option_type == "call":
return S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
elif option_type == "put":
return K * np.exp(-r * T) * norm.cdf(-d2) - S * norm.cdf(-d1)
else:
raise ValueError("Invalid option type. Use 'call' or 'put'.")
# Greeks calculation functions
def calculate_greeks(S, K, T, r, sigma, option_type="call"):
d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
d2 = d1 - sigma * np.sqrt(T)
delta = norm.cdf(d1) if option_type == "call" else -norm.cdf(-d1)
gamma = norm.pdf(d1) / (S * sigma * np.sqrt(T))
vega = S * norm.pdf(d1) * np.sqrt(T) / 100 # Expressed per 1% change in volatility
theta = - (S * norm.pdf(d1) * sigma) / (2 * np.sqrt(T))
if option_type == "call":
theta -= r * K * np.exp(-r * T) * norm.cdf(d2)
else:
theta += r * K * np.exp(-r * T) * norm.cdf(-d2)
theta /= 365 # Per day
rho = K * T * np.exp(-r * T) * (norm.cdf(d2) if option_type == "call" else -norm.cdf(-d2)) / 100 # Per 1% change in rates
return delta, gamma, vega, theta, rho
# Load the Excel file
file_path = r'PUT YOUR PERSONAL DIRECTORY'
# Load the second sheet of the Excel file
data = pd.ExcelFile(file_path)
df = data.parse(data.sheet_names[1]) # Load the second sheet
# Data cleaning
columns_mapping = {
"Unnamed: 0": "Date",
"ISP IM Equity": "ISP_Price",
"UCG IM Equity": "UCG_Price",
"SAN SM Equity": "SAN_Price",
"BBVA SM Equity": "BBVA_Price",
"ACA FP Equity": "ACA_Price",
"BNP FP Equity": "BNP_Price",
"EUDR1T Curncy": "RiskFreeRate"
}
df = df.rename(columns=columns_mapping)
df = df[~df["Date"].astype(str).str.contains("DATES", na=False)]
df["Date"] = pd.to_datetime(df["Date"], errors='coerce')
df = df.dropna(subset=["Date", "RiskFreeRate"])
df.set_index("Date", inplace=True)
stocks = ["ISP_Price", "UCG_Price", "SAN_Price", "BBVA_Price", "ACA_Price", "BNP_Price"]
for stock in stocks:
df_stock = df[[stock, "RiskFreeRate"]].dropna()
quarterly_strike_prices = df_stock[stock].resample('Q').mean()
print(f"\nQuarterly Average Prices (Strike Prices) for {stock}:")
print(quarterly_strike_prices)
# Calculation of parameters for Black-Scholes
current_price = df_stock[stock].iloc[-1]
strike_price = quarterly_strike_prices.iloc[-1]
time_to_maturity = 0.25 # Quarter in years
#We assumed a strike price #given by the quarterly moving average price of the underlying asset
risk_free_rate = df_stock["RiskFreeRate"].iloc[-1] / 100 # Convert to decimal form
volatility = df_stock[stock].pct_change().std() * np.sqrt(252) # Annualized volatility
# Calculation of option prices
call_price = black_scholes(current_price, strike_price, time_to_maturity, risk_free_rate, volatility, option_type="call")
put_price = black_scholes(current_price, strike_price, time_to_maturity, risk_free_rate, volatility, option_type="put")
# Output of results
print(f"Call Option Price for {stock}: {call_price:.2f}")
print(f"Put Option Price for {stock}: {put_price:.2f}")
# Creating the price chart and quarterly strike prices
plt.figure(figsize=(10, 6))
# Plot the historical stock price (blue line)
plt.plot(df_stock.index, df_stock[stock], label=f'Price {stock}', color='blue')
# Plot the Call and Put option prices
plt.plot(df_stock.index, [black_scholes(price, strike_price, time_to_maturity, risk_free_rate, volatility, option_type="call") for price in df_stock[stock]],
label=f'Call Price {stock}', color='green', linestyle='--')
plt.plot(df_stock.index, [black_scholes(price, strike_price, time_to_maturity, risk_free_rate, volatility, option_type="put") for price in df_stock[stock]],
label=f'Put Price {stock}', color='red', linestyle='--')
# Plot the strike price (orange dashed line)
plt.plot(quarterly_strike_prices.index, quarterly_strike_prices, label='Strike Price (Quarterly Average)', color='orange', linestyle='--')
# Add title, labels, and legend
plt.title(f'Trend of {stock} Price and Option Prices')
plt.xlabel('Date')
plt.ylabel('Price')
plt.legend()
plt.grid(True)
# Show the chart
plt.show()
Quarterly Average Prices (Strike Prices) for ISP_Price: Date 2022-03-31 2.399073 2022-06-30 1.944486 2022-09-30 1.756502 2022-12-31 1.997592 2023-03-31 2.380092 2023-06-30 2.366968 2023-09-30 2.462891 2023-12-31 2.542603 Freq: Q-DEC, Name: ISP_Price, dtype: float64 Call Option Price for ISP_Price: 0.24 Put Option Price for ISP_Price: 0.11
Quarterly Average Prices (Strike Prices) for UCG_Price: Date 2022-03-31 12.527375 2022-06-30 9.711603 2022-09-30 9.712492 2022-12-31 12.398500 2023-03-31 17.019292 2023-06-30 18.972726 2023-09-30 22.032812 2023-12-31 23.888016 Freq: Q-DEC, Name: UCG_Price, dtype: float64 Call Option Price for UCG_Price: 2.48 Put Option Price for UCG_Price: 1.57
Quarterly Average Prices (Strike Prices) for SAN_Price: Date 2022-03-31 3.116875 2022-06-30 2.881230 2022-09-30 2.508023 2022-12-31 2.682164 2023-03-31 3.342608 2023-06-30 3.261331 2023-09-30 3.515046 2023-12-31 3.672008 Freq: Q-DEC, Name: SAN_Price, dtype: float64 Call Option Price for SAN_Price: 0.32 Put Option Price for SAN_Price: 0.17
Quarterly Average Prices (Strike Prices) for BBVA_Price: Date 2022-03-31 5.465016 2022-06-30 4.764548 2022-09-30 4.527530 2022-12-31 5.297406 2023-03-31 6.629092 2023-06-30 6.565645 2023-09-30 7.150738 2023-12-31 8.021778 Freq: Q-DEC, Name: BBVA_Price, dtype: float64 Call Option Price for BBVA_Price: 0.66 Put Option Price for BBVA_Price: 0.38
Quarterly Average Prices (Strike Prices) for ACA_Price: Date 2022-03-31 12.328250 2022-06-30 9.991905 2022-09-30 9.096636 2022-12-31 9.319141 2023-03-31 10.857677 2023-06-30 11.034161 2023-09-30 11.334369 2023-12-31 11.891270 Freq: Q-DEC, Name: ACA_Price, dtype: float64 Call Option Price for ACA_Price: 1.36 Put Option Price for ACA_Price: 0.28
Quarterly Average Prices (Strike Prices) for BNP_Price: Date 2022-03-31 58.730234 2022-06-30 50.353333 2022-09-30 46.458182 2022-12-31 50.005313 2023-03-31 59.996615 2023-06-30 57.141774 2023-09-30 58.876769 2023-12-31 58.287302 Freq: Q-DEC, Name: BNP_Price, dtype: float64 Call Option Price for BNP_Price: 6.66 Put Option Price for BNP_Price: 1.78
import pandas as pd
import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt
from scipy.interpolate import make_interp_spline
# Black-Scholes function to calculate option price
def black_scholes(S, K, T, r, sigma, option_type="call"):
d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
d2 = d1 - sigma * np.sqrt(T)
if option_type == "call":
return S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
elif option_type == "put":
return K * np.exp(-r * T) * norm.cdf(-d2) - S * norm.cdf(-d1)
else:
raise ValueError("Invalid option type. Use 'call' or 'put'.")
# Function to calculate the Greeks
def calculate_greeks(S, K, T, r, sigma, option_type="call"):
d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
d2 = d1 - sigma * np.sqrt(T)
delta = norm.cdf(d1) if option_type == "call" else -norm.cdf(-d1)
gamma = norm.pdf(d1) / (S * sigma * np.sqrt(T))
vega = S * norm.pdf(d1) * np.sqrt(T) / 100 # Vega per 1% change in volatility
theta = - (S * norm.pdf(d1) * sigma) / (2 * np.sqrt(T))
if option_type == "call":
theta -= r * K * np.exp(-r * T) * norm.cdf(d2)
else:
theta += r * K * np.exp(-r * T) * norm.cdf(-d2)
theta /= 365 # Convert theta to per-day value
rho = K * T * np.exp(-r * T) * (norm.cdf(d2) if option_type == "call" else -norm.cdf(-d2)) / 100 # Rho per 1% change in rates
return delta, gamma, vega, theta, rho
# Load the Excel file
file_path = r'PUT YOUR PERSONAL DIRECTORY'
# Load the second sheet from the Excel file
data = pd.ExcelFile(file_path)
df = data.parse(data.sheet_names[1]) # Load second sheet
# Data cleaning and renaming columns
columns_mapping = {
"Unnamed: 0": "Date",
"ISP IM Equity": "ISP_Price",
"UCG IM Equity": "UCG_Price",
"SAN SM Equity": "SAN_Price",
"BBVA SM Equity": "BBVA_Price",
"ACA FP Equity": "ACA_Price",
"BNP FP Equity": "BNP_Price",
"EUDR1T Curncy": "RiskFreeRate"
}
df = df.rename(columns=columns_mapping)
df = df[~df["Date"].astype(str).str.contains("DATES", na=False)]
df["Date"] = pd.to_datetime(df["Date"], errors='coerce')
df = df.dropna(subset=["Date", "RiskFreeRate"])
df.set_index("Date", inplace=True)
# List of stocks to analyze
stocks = ["ISP_Price", "UCG_Price", "SAN_Price", "BBVA_Price", "ACA_Price", "BNP_Price"]
for stock in stocks:
df_stock = df[[stock, "RiskFreeRate"]].dropna()
quarterly_strike_prices = df_stock[stock].resample('Q').mean()
# Compute parameters for Black-Scholes model
current_price = df_stock[stock].iloc[-1]
strike_price = quarterly_strike_prices.iloc[-1]
time_to_maturity = 0.25 # One quarter (3 months) in years
risk_free_rate = df_stock["RiskFreeRate"].iloc[-1] / 100 # Convert percentage to decimal
volatility = df_stock[stock].pct_change().std() * np.sqrt(252) # Annualized volatility
# Generate a range of underlying asset prices
prices = np.linspace(current_price * 0.7, current_price * 1.3, 100) # Range from 70% to 130% of the current price
# Calculate Greeks for each price in the range
deltas, gammas, vegas, thetas, rhos = [], [], [], [], []
for price in prices:
delta, gamma, vega, theta, rho = calculate_greeks(price, strike_price, time_to_maturity, risk_free_rate, volatility, option_type="call")
deltas.append(delta)
gammas.append(gamma)
vegas.append(vega)
thetas.append(theta)
rhos.append(rho)
# Plot Greeks against stock price
plt.figure(figsize=(15, 10))
# Delta vs Price
plt.subplot(2, 3, 1)
plt.plot(prices, deltas, label="Delta", color='blue')
plt.title(f"Delta vs Price for {stock}")
plt.xlabel('Underlying Price')
plt.ylabel('Delta')
plt.grid(True)
# Gamma vs Price
plt.subplot(2, 3, 2)
plt.plot(prices, gammas, label="Gamma", color='green')
plt.title(f"Gamma vs Price for {stock}")
plt.xlabel('Underlying Price')
plt.ylabel('Gamma')
plt.grid(True)
# Vega vs Price
plt.subplot(2, 3, 3)
plt.plot(prices, vegas, label="Vega", color='purple')
plt.title(f"Vega vs Price for {stock}")
plt.xlabel('Underlying Price')
plt.ylabel('Vega')
plt.grid(True)
# Theta vs Price
plt.subplot(2, 3, 4)
plt.plot(prices, thetas, label="Theta", color='red')
plt.title(f"Theta vs Price for {stock}")
plt.xlabel('Underlying Price')
plt.ylabel('Theta')
plt.grid(True)
# Rho vs Price
plt.subplot(2, 3, 5)
plt.plot(prices, rhos, label="Rho", color='orange')
plt.title(f"Rho vs Price for {stock}")
plt.xlabel('Underlying Price')
plt.ylabel('Rho')
plt.grid(True)
plt.tight_layout()
plt.show()
import pandas as pd
import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt
# Black-Scholes pricing function
def black_scholes(S, K, T, r, sigma, option_type="call"):
d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
d2 = d1 - sigma * np.sqrt(T)
if option_type == "call":
return S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
elif option_type == "put":
return K * np.exp(-r * T) * norm.cdf(-d2) - S * norm.cdf(-d1)
else:
raise ValueError("Invalid option type. Use 'call' or 'put'.")
# Assuming 'df' is already cleaned and contains necessary columns
stocks = ["ISP_Price", "UCG_Price", "SAN_Price", "BBVA_Price", "ACA_Price", "BNP_Price"]
for stock in stocks:
df_stock = df[[stock, "RiskFreeRate"]].dropna()
quarterly_strike_prices = df_stock[stock].resample('Q').mean()
# Extract parameters
S = df_stock[stock].iloc[-1] # Current price
r = df_stock["RiskFreeRate"].iloc[-1] / 100 # Risk-free rate (decimal)
sigma = df_stock[stock].pct_change().std() * np.sqrt(252) # Annualized volatility
# Generate option prices vs. strike price (K)
K_values = np.linspace(S * 0.8, S * 1.2, 50)
option_prices_vs_K_call = [black_scholes(S, K, 1, r, sigma, option_type="call") for K in K_values]
option_prices_vs_K_put = [black_scholes(S, K, 1, r, sigma, option_type="put") for K in K_values]
# Generate option prices vs. time to maturity (T)
T_values = np.linspace(0.01, 2, 50)
option_prices_vs_T_call = [black_scholes(S, S, T, r, sigma, option_type="call") for T in T_values]
option_prices_vs_T_put = [black_scholes(S, S, T, r, sigma, option_type="put") for T in T_values]
# Create a figure with subplots for better layout
fig, ax = plt.subplots(2, 1, figsize=(12, 12))
# Plot option prices vs. strike price (K)
ax[0].plot(K_values, option_prices_vs_K_call, label=f'{stock} Call Price vs Strike Price', color='blue', linestyle='-')
ax[0].plot(K_values, option_prices_vs_K_put, label=f'{stock} Put Price vs Strike Price', color='red', linestyle='--')
ax[0].set_title(f'{stock} - Option Price vs Strike Price', fontsize=14)
ax[0].set_xlabel('Strike Price (K)', fontsize=12)
ax[0].set_ylabel('Option Price', fontsize=12)
ax[0].grid(True)
ax[0].legend(loc='best')
# Plot option prices vs. time to maturity (T)
ax[1].plot(T_values, option_prices_vs_T_call, label=f'{stock} Call Price vs Time to Maturity', color='green', linestyle='-')
ax[1].plot(T_values, option_prices_vs_T_put, label=f'{stock} Put Price vs Time to Maturity', color='orange', linestyle='--')
ax[1].set_title(f'{stock} - Option Price vs Time to Maturity', fontsize=14)
ax[1].set_xlabel('Time to Maturity (T)', fontsize=12)
ax[1].set_ylabel('Option Price', fontsize=12)
ax[1].grid(True)
ax[1].legend(loc='best')
# Adjust the layout to prevent overlap and display the plots
plt.tight_layout()
plt.show()
import pandas as pd
import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# Function to calculate implied volatility
def implied_volatility(option_price, S, K, T, r, option_type="call", tol=1e-5, max_iter=100):
sigma = 0.2 # Initial guess for volatility
for i in range(max_iter):
price = black_scholes(S, K, T, r, sigma, option_type)
vega = S * norm.pdf((np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))) * np.sqrt(T)
if vega < 1e-6: # Avoid division by zero
return np.nan
price_diff = option_price - price
if abs(price_diff) < tol:
return sigma
sigma += price_diff / vega
return np.nan # Return NaN if no convergence
# Black-Scholes pricing function
def black_scholes(S, K, T, r, sigma, option_type="call"):
if T <= 0 or sigma <= 0:
return 0.0 # Return 0 for invalid inputs
d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
d2 = d1 - sigma * np.sqrt(T)
if option_type == "call":
return S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
elif option_type == "put":
return K * np.exp(-r * T) * norm.cdf(-d2) - S * norm.cdf(-d1)
else:
raise ValueError("Invalid option type. Use 'call' or 'put'.")
# Load the data from the existing dataframe
# Assuming 'df' contains the cleaned and structured data as in the previous code
stocks = ["ISP_Price", "UCG_Price", "SAN_Price", "BBVA_Price", "ACA_Price", "BNP_Price"]
for stock in stocks:
df_stock = df[[stock, "RiskFreeRate"]].dropna()
quarterly_strike_prices = df_stock[stock].resample('Q').mean()
print(f"\nCalculating implied volatility surface for {stock}:")
# Extract parameters
S = df_stock[stock].iloc[-1] # Current price
r = df_stock["RiskFreeRate"].iloc[-1] / 100 # Risk-free rate (decimal)
# Compute historical volatility (annualized)
daily_returns = df_stock[stock].pct_change().dropna()
historical_volatility = daily_returns.std() * np.sqrt(252) # Annualized volatility
# Debugging step: print key parameters
print(f"Stock: {stock}, Current Price: {S}, Risk-free Rate: {r:.4f}, Historical Volatility: {historical_volatility:.4f}")
print("Quarterly Strike Prices:")
print(quarterly_strike_prices)
# Grid for strike prices and time to maturity
K_values = np.linspace(S * 0.8, S * 1.2, 30)
T_values = np.linspace(0.1, 2, 30)
K_grid, T_grid = np.meshgrid(K_values, T_values)
# Compute implied volatilities for the surface
IV_surface = np.zeros_like(K_grid)
for i in range(K_grid.shape[0]):
for j in range(K_grid.shape[1]):
market_price = black_scholes(S, K_grid[i, j], T_grid[i, j], r, historical_volatility, option_type="call")
IV_surface[i, j] = implied_volatility(market_price, S, K_grid[i, j], T_grid[i, j], r, option_type="call")
# Replace NaN values with a small constant for visualization
IV_surface = np.nan_to_num(IV_surface, nan=0.01)
# Plot implied volatility surface
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(K_grid, T_grid, IV_surface, cmap='viridis')
ax.set_title(f"Implied Volatility Surface for {stock}")
ax.set_xlabel("Strike Price (K)")
ax.set_ylabel("Time to Maturity (T)")
ax.set_zlabel("Implied Volatility")
plt.show()
Calculating implied volatility surface for ISP_Price: Stock: ISP_Price, Current Price: 2.6435, Risk-free Rate: 0.0398, Historical Volatility: 0.3212 Quarterly Strike Prices: Date 2022-03-31 2.399073 2022-06-30 1.944486 2022-09-30 1.756502 2022-12-31 1.997592 2023-03-31 2.380092 2023-06-30 2.366968 2023-09-30 2.462891 2023-12-31 2.542603 Freq: Q-DEC, Name: ISP_Price, dtype: float64
Calculating implied volatility surface for UCG_Price: Stock: UCG_Price, Current Price: 24.565, Risk-free Rate: 0.0398, Historical Volatility: 0.4156 Quarterly Strike Prices: Date 2022-03-31 12.527375 2022-06-30 9.711603 2022-09-30 9.712492 2022-12-31 12.398500 2023-03-31 17.019292 2023-06-30 18.972726 2023-09-30 22.032812 2023-12-31 23.888016 Freq: Q-DEC, Name: UCG_Price, dtype: float64
Calculating implied volatility surface for SAN_Price: Stock: SAN_Price, Current Price: 3.7795, Risk-free Rate: 0.0398, Historical Volatility: 0.3215 Quarterly Strike Prices: Date 2022-03-31 3.116875 2022-06-30 2.881230 2022-09-30 2.508023 2022-12-31 2.682164 2023-03-31 3.342608 2023-06-30 3.261331 2023-09-30 3.515046 2023-12-31 3.672008 Freq: Q-DEC, Name: SAN_Price, dtype: float64
Calculating implied volatility surface for BBVA_Price: Stock: BBVA_Price, Current Price: 8.226, Risk-free Rate: 0.0398, Historical Volatility: 0.3158 Quarterly Strike Prices: Date 2022-03-31 5.465016 2022-06-30 4.764548 2022-09-30 4.527530 2022-12-31 5.297406 2023-03-31 6.629092 2023-06-30 6.565645 2023-09-30 7.150738 2023-12-31 8.021778 Freq: Q-DEC, Name: BBVA_Price, dtype: float64
Calculating implied volatility surface for ACA_Price: Stock: ACA_Price, Current Price: 12.852, Risk-free Rate: 0.0398, Historical Volatility: 0.2803 Quarterly Strike Prices: Date 2022-03-31 12.328250 2022-06-30 9.991905 2022-09-30 9.096636 2022-12-31 9.319141 2023-03-31 10.857677 2023-06-30 11.034161 2023-09-30 11.334369 2023-12-31 11.891270 Freq: Q-DEC, Name: ACA_Price, dtype: float64
Calculating implied volatility surface for BNP_Price: Stock: BNP_Price, Current Price: 62.59, Risk-free Rate: 0.0398, Historical Volatility: 0.3107 Quarterly Strike Prices: Date 2022-03-31 58.730234 2022-06-30 50.353333 2022-09-30 46.458182 2022-12-31 50.005313 2023-03-31 59.996615 2023-06-30 57.141774 2023-09-30 58.876769 2023-12-31 58.287302 Freq: Q-DEC, Name: BNP_Price, dtype: float64
import pandas as pd
import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt
from scipy.interpolate import griddata
# Black-Scholes pricing function
def black_scholes(S, K, T, r, sigma, option_type="call"):
d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
d2 = d1 - sigma * np.sqrt(T)
if option_type == "call":
return S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
elif option_type == "put":
return K * np.exp(-r * T) * norm.cdf(-d2) - S * norm.cdf(-d1)
else:
raise ValueError("Invalid option type. Use 'call' or 'put'.")
# Function to calculate the local volatility using Dupire's formula (approximation)
def calculate_local_volatility(S, K, T, r, sigma, epsilon=1e-5):
# Calculate the Black-Scholes price for the option
C = black_scholes(S, K, T, r, sigma, option_type="call")
# Calculate the derivative of the option price with respect to time to maturity (T)
C_T_plus = black_scholes(S, K, T + epsilon, r, sigma, option_type="call")
C_T_minus = black_scholes(S, K, T - epsilon, r, sigma, option_type="call")
dC_dT = (C_T_plus - C_T_minus) / (2 * epsilon) # Numerical approximation
# Calculate the second derivative of the option price with respect to the stock price (S)
C_S_plus = black_scholes(S + epsilon, K, T, r, sigma, option_type="call")
C_S_minus = black_scholes(S - epsilon, K, T, r, sigma, option_type="call")
d2C_dS2 = (C_S_plus - 2 * C + C_S_minus) / (epsilon ** 2) # Numerical approximation
# Local volatility using the Dupire formula approximation
local_vol = np.sqrt(2 * dC_dT / (S ** 2 * d2C_dS2))
return local_vol
# Load the data from the existing dataframe
# Assuming 'df' contains the cleaned and structured data as in the previous code
stocks = ["ISP_Price", "UCG_Price", "SAN_Price", "BBVA_Price", "ACA_Price", "BNP_Price"]
for stock in stocks:
df_stock = df[[stock, "RiskFreeRate"]].dropna()
print(f"\nCalculating Local Volatility Surface for {stock}:")
# Extract parameters
S_current = df_stock[stock].iloc[-1] # Current price
r = df_stock["RiskFreeRate"].iloc[-1] / 100 # Risk-free rate
sigma = df_stock[stock].pct_change().std() * np.sqrt(252) # Volatility
# Debugging step: print key parameters
print(f"Stock: {stock}, Current Price: {S_current}, Risk-free Rate: {r:.4f}, Volatility: {sigma:.4f}")
# Grid for stock prices and time to maturity
K_values = np.linspace(S_current * 0.8, S_current * 1.2, 30)
T_values = np.linspace(0.1, 2, 30)
K_grid, T_grid = np.meshgrid(K_values, T_values)
# Compute local volatilities for the surface
local_vol_surface = np.zeros_like(K_grid)
for i in range(K_grid.shape[0]):
for j in range(K_grid.shape[1]):
local_vol_surface[i, j] = calculate_local_volatility(S_current, K_grid[i, j], T_grid[i, j], r, sigma)
# Plot local volatility surface
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(K_grid, T_grid, local_vol_surface, cmap='viridis')
ax.set_title(f"Local Volatility Surface for {stock}")
ax.set_xlabel("Strike Price (K)")
ax.set_ylabel("Time to Maturity (T)")
ax.set_zlabel("Local Volatility")
plt.show()
Calculating Local Volatility Surface for ISP_Price: Stock: ISP_Price, Current Price: 2.6435, Risk-free Rate: 0.0398, Volatility: 0.3212
Calculating Local Volatility Surface for UCG_Price: Stock: UCG_Price, Current Price: 24.565, Risk-free Rate: 0.0398, Volatility: 0.4156
Calculating Local Volatility Surface for SAN_Price: Stock: SAN_Price, Current Price: 3.7795, Risk-free Rate: 0.0398, Volatility: 0.3215
Calculating Local Volatility Surface for BBVA_Price: Stock: BBVA_Price, Current Price: 8.226, Risk-free Rate: 0.0398, Volatility: 0.3158
Calculating Local Volatility Surface for ACA_Price: Stock: ACA_Price, Current Price: 12.852, Risk-free Rate: 0.0398, Volatility: 0.2803
Calculating Local Volatility Surface for BNP_Price: Stock: BNP_Price, Current Price: 62.59, Risk-free Rate: 0.0398, Volatility: 0.3107
import pandas as pd
import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# Black-Scholes pricing function
def black_scholes(S, K, T, r, sigma, option_type="call"):
d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
d2 = d1 - sigma * np.sqrt(T)
if option_type == "call":
return S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
elif option_type == "put":
return K * np.exp(-r * T) * norm.cdf(-d2) - S * norm.cdf(-d1)
else:
raise ValueError("Invalid option type. Use 'call' or 'put'.")
# Function to calculate local volatility using Dupire's formula (numerically approximated)
def calculate_local_volatility(S, K, T, r, sigma, epsilon=1e-5):
# Calculate the Black-Scholes price for the option
C = black_scholes(S, K, T, r, sigma, option_type="call")
# Calculate the derivative of the option price with respect to time to maturity (T)
C_T_plus = black_scholes(S, K, T + epsilon, r, sigma, option_type="call")
C_T_minus = black_scholes(S, K, T - epsilon, r, sigma, option_type="call")
dC_dT = (C_T_plus - C_T_minus) / (2 * epsilon) # Numerical approximation
# Calculate the second derivative of the option price with respect to the stock price (S)
C_S_plus = black_scholes(S + epsilon, K, T, r, sigma, option_type="call")
C_S_minus = black_scholes(S - epsilon, K, T, r, sigma, option_type="call")
d2C_dS2 = (C_S_plus - 2 * C + C_S_minus) / (epsilon ** 2) # Numerical approximation
# Local volatility using the Dupire formula approximation
local_vol = np.sqrt(2 * dC_dT / (S ** 2 * d2C_dS2))
return local_vol
# Function to calculate implied volatility using Newton-Raphson method
def implied_volatility(option_price, S, K, T, r, option_type="call", tol=1e-5, max_iter=100):
sigma = 0.2 # Initial guess for volatility
for _ in range(max_iter):
price = black_scholes(S, K, T, r, sigma, option_type)
vega = S * norm.pdf((np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))) * np.sqrt(T)
if vega < 1e-6:
return np.nan # Avoid division by zero
price_diff = option_price - price
if abs(price_diff) < tol:
return sigma
sigma += price_diff / vega
return np.nan # Return NaN if no convergence
# Load the data from the existing dataframe
# Assuming 'df' contains the cleaned and structured data as in the previous code
stocks = ["ISP_Price", "UCG_Price", "SAN_Price", "BBVA_Price", "ACA_Price", "BNP_Price"]
for stock in stocks:
df_stock = df[[stock, "RiskFreeRate"]].dropna()
print(f"\nCalculating Local vs Implied Volatility for {stock}:")
# Extract parameters
S_current = df_stock[stock].iloc[-1] # Current price
r = df_stock["RiskFreeRate"].iloc[-1] / 100 # Risk-free rate
sigma = df_stock[stock].pct_change().std() * np.sqrt(252) # Volatility
# Debugging step: print key parameters
print(f"Stock: {stock}, Current Price: {S_current}, Risk-free Rate: {r:.4f}, Volatility: {sigma:.4f}")
# Grid for stock prices and time to maturity
K_values = np.linspace(S_current * 0.8, S_current * 1.2, 30)
T_values = np.linspace(0.1, 2, 30)
K_grid, T_grid = np.meshgrid(K_values, T_values)
# Compute local volatilities and implied volatilities for the surface
local_vol_surface = np.zeros_like(K_grid)
implied_vol_surface = np.zeros_like(K_grid)
for i in range(K_grid.shape[0]):
for j in range(K_grid.shape[1]):
market_price = black_scholes(S_current, K_grid[i, j], T_grid[i, j], r, sigma, option_type="call")
implied_vol_surface[i, j] = implied_volatility(market_price, S_current, K_grid[i, j], T_grid[i, j], r, option_type="call")
local_vol_surface[i, j] = calculate_local_volatility(S_current, K_grid[i, j], T_grid[i, j], r, sigma)
# Replace NaN values in the implied volatility surface
implied_vol_surface = np.nan_to_num(implied_vol_surface, nan=0.01)
# Plot implied volatility surface
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(121, projection='3d')
ax.plot_surface(K_grid, T_grid, implied_vol_surface, cmap='viridis')
ax.set_title(f"Implied Volatility Surface for {stock}")
ax.set_xlabel("Strike Price (K)")
ax.set_ylabel("Time to Maturity (T)")
ax.set_zlabel("Volatility")
# Plot local volatility surface
ax2 = fig.add_subplot(122, projection='3d')
ax2.plot_surface(K_grid, T_grid, local_vol_surface, cmap='plasma')
ax2.set_title(f"Local Volatility Surface for {stock}")
ax2.set_xlabel("Strike Price (K)")
ax2.set_ylabel("Time to Maturity (T)")
ax2.set_zlabel("Volatility")
plt.tight_layout()
plt.show()
Calculating Local vs Implied Volatility for ISP_Price: Stock: ISP_Price, Current Price: 2.6435, Risk-free Rate: 0.0398, Volatility: 0.3212
Calculating Local vs Implied Volatility for UCG_Price: Stock: UCG_Price, Current Price: 24.565, Risk-free Rate: 0.0398, Volatility: 0.4156
Calculating Local vs Implied Volatility for SAN_Price: Stock: SAN_Price, Current Price: 3.7795, Risk-free Rate: 0.0398, Volatility: 0.3215
Calculating Local vs Implied Volatility for BBVA_Price: Stock: BBVA_Price, Current Price: 8.226, Risk-free Rate: 0.0398, Volatility: 0.3158
Calculating Local vs Implied Volatility for ACA_Price: Stock: ACA_Price, Current Price: 12.852, Risk-free Rate: 0.0398, Volatility: 0.2803
Calculating Local vs Implied Volatility for BNP_Price: Stock: BNP_Price, Current Price: 62.59, Risk-free Rate: 0.0398, Volatility: 0.3107